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PSEUDOSPECTRAL  SOLUTION  OF  ONE  DIMENSIONAL  AND  TWO  DIMENSIONAL 
INVISCID  FLOWS  WITH  SHOCK  WAVES 


I.  INTRODUCTION 

Pseud ospectral  techniques  have  been  used  to  solve  the  one  dimensional 
propagating  shock  wave  problem.  Taylor  et  al  (Reference  I)  and  Gottlieb  et 
al  (Reference  2)  have  done  so  using  the  Euler  equations  of  motion.  Taylor 
utilized  the  FCT  (Flux  Corrected  Transport)  algorithm  of  Boris  and  Book 
(Reference  3)  to  damp  out  unwanted  numerical  oscillations.  This  procedure 
yielded  broadening  of  the  shock  wave.  They  treated  a  Mach  1.4  shock  wave 
propagating  into  a  free  stream  at  rest.  The  flow  behind  the  shock  wave  was 
subsonic.  Gottlieb  et  al  treated  the  shock  tube  problem  for  shock  wave  Mach 
nunbers  of  2.1  and  29.3.  The  free  stream  was  subsonic  with  the  flow  behind 
the  shock  wave  being  supersonic  for  both  roach  number  cases.  Ihey  performed 
a  detailed  analysis  of  the  effects  of  different  filtering  techniques  on 
reducing  unwanted  numerical  oscillations.  They  considered  the  Shuman  filter 
given  by: 

ri  ■  ”3  +  V*  &  3  1 +  V  "5-i 1  <l) 

is  the  filtered  conservative  variable  at  the  spatial  location  and  the 
nctl  time  step.  For  a  two  dimensional  problem  one  would  need  to  filter  in 
each  direction  separately.  The  0.J,  coefficients  are  given  by 
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e  .  .  giSttLl  fatililSttLl  pj),t.(pj'  pj-.Qi 
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and  8  is  a  constant  greater  than  zero  and  less  than  one*  Beta  was  chosen  to 
be  0.01.  The  above  was  used  in  one  of  two  versions, 

constant  0^  coefficients  and  variable  0^  ^coefficients*  The  former  is 
qualitatively  equivalent  to  a  first  order  artificial  viscosity  scheme*  Both 
were  applied  to  the  physical  variables  directly*  They  also  utilized  a  low 
pass  spectral  filter,  which  they  developed,  to  damp  out  the  oscillations 
which  arose  from  the  highest  frequency  spectral  components*  The  form  of 
their  spectral  filter  is: 


e 


k-k 


°fc- 


max 


(3) 


where  k.  is  the  spectral  wavenumber,  kmax  is  the  maxinum  wavenumber 
corresponding  to  the  total  number  of  collocation  points  and  kQ*^  N  where  N 
is  the  total  number  of  collocation  points  used  to  represent  the  flow*  The 
spectral  filter  was  applied  first,  followed  by  the  Shuman  filter*  They 
determined  rules  for  applying  the  low  pass  spectral  filter.  They  found  that 
applying  it  over  the  highest  sixth  of  the  frequency  values  gave  good 
results.  The  Shuman  filtering  employed  was  one  sided.  That  is,  the  shock 
position  was  determined  first  and  then  the  filter  was  applied  over  the 
region  behind  the  shock  wave  and  separately  to  the  region  in  front  of  the 
shock  wave*  Using  this  approach  they  were  able  to  obtain  a  sharp  shock  with 
the  correct  propagation  velocity.  Both  approaches,  however,  have  some 
dratbadcs.  The  former  did  not  yield  a  sharp  discontinuity  while  the  latter 
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required  an  examination  of  the  spectral  coefficients  at  each  time  step  to 
determine  the  shock  wave  location  in  order  to  avoid  applying  the  physical 
space  filter  across  the  shock  front.  When  more  general  classes  of  inviscid 
flows  are  treated  (ones  with  complex,  multiple  shock  geometries)  the 
smearing  in  the  first  approach  may  prove  unacceptable.  The  one-sided 
smoothing  of  the  latter  will  become  cunbersome  to  employ. 

A  brief  outline  of  pseud  os  pectral  techniques  will  be  given  in 
Section  2 •  The  third  section  of  this  report  will  present  results  for  the 
one  dimensional  propagating  shodc  wave  problem  using  a  different  physical 
space  smoothing  function  than  either  oftjj while  retaining  the 
lowpass  spectral  f iltering  t^ahrtflCque  of  Reference  2.  M  artificial 
viscosity  scheme  is  used  uniformly  throughout  the  entire  flow  field, 
including  across  the  shock  front,  to  resolve  the  shock  wave  as  a  sharp 
discontinuity  and  at  the  same  time  maintain  the  correct  shock  propagation 
velocity. 

To  further  demonstrate  the  utility  of  this  approach  to  the  solution  of 
flows  by  pseudospectral  methods,  solutions  to  two-dimensional  inviscid 
supersonic  wedge  flows  will  also  be  presented  in  Section  4  of  this  report. 

To  the  present  author's  knowledge,  this  is  the  first  time  pseudospectral 
solution  techniques  have  been  used  to  successfully  treat  two-dimensional 
inviscid  flows. 

2.  PSEUDOSPECTRAL  METHODS 

A  brief  description  of  pseud  os pectral  techniques  will  be  presented  here 
for  completeness.  For  those  readers  Interested  in  a  detailed  exposition  on 
pseudospectral  techniques.  Reference  4  is  strongly  recommended. 
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Pseudospectral  techniques  involve  the  use  of  series  of  functions  to 
solve  differential  equations.  For  all  work  reported  herein,  Che by she v 
polynomials  are  used.  Chetyshev  polynomials  are  represented  by  TQ(x)  and 
are  given  by: 

Tn(x)  ■  cos  [n  cos“*(x) ] 

-  cos  [n  9  ]  (4a) 

where  9  *  cos”^  (x) 


The  Chetyshev  polynomials  have  the  following  property: 


2  TnU) 


n+1 


T'n-l(x) 

n-1 


(4b) 


Chebyshev  polynomials  may  be  used  to  represent  a  function  F(x)  in  the 
following  manner 


N 

F(x)  -  I  AnTn(x)  (5) 

n-o 


A  function  F  (x,t)  may  be  represented  as: 


N 

F  (x *  t )  -  Z  An(t)  T  (x)  (6) 

n-o 


where  the  time  dependence  is  totally  contained  in  the  series  coefficients 
An(t),  and  the  x  dependence  in  the  Chebyshev  functions  Tn(x). 
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Let  us  now  consider  the  application  of  such  techniques  for  the  solution 
of  the  onetime ns ional  unsteady  Euler  equation  in  conservation  law  form: 

where 

and 


The  pseudospectral  solution,  using  Chetyshev  collocation,  of  this  equation 
involves  using  Chebyshev  series  to  obtain  the  spatial  derivative  and  finite 
difference  algorithms  to  obtain  the  time  derivative.  (A  flowchart  of  the 
solution  procedure  is  shown  in  Figure  1).  Collocation  involves  the 
specification  of  the  initial  flow  variables  and  the  conputation  of  the  time 
dependent  solution  for  the  flow  variables  at  distinct  pre -determined  spatial 
positions  or  points.  These  positions  are  the  collocation  points.  The 
spatial  derivative  is  obtained  as  follows.  At  tQ  the  values  of  eT(x)  at  the 
collocation  points  Xj  are  specified. 

The  collocation  points  are  given  by 

Xj-  cos  ySL)  0<  j  <  N  (8) 

where  N  is  the  total  number  of  Chetyshev  polynomials  one  chooses  to  use  to 
represent  the  function  eT(x)  •  As  can  easily  be  seen,  the  Xj  points  are  not 
evenly  spaced,  but  are  clustered  about  x  ■  ±  1. 

We  represent  lf(x)  by 

N 

f(t,x)  -  I  An(t)Tn(x)  (9) 


JU  J?.  o 

at  fc 


(7  a) 


p 

2  " 

pi 

E  - 

p+pu 

e 

(e+p)u 

P 

,  2 

(r  1) 


(7b) 


(7c) 
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The  left  hand  side  of  this  equation  is  known  while  the  are  at  this 

point  unknown.  The  first  step  is  therefore  to  solve  for  the  A^s  for  each 
£(x)  vector  element.  This  could  be  done  by  a  simple  matrix  inversion. 
However,  it  is  much  faster  to  use  FFT's  (fast  Fourier  transforms).  Vfe 
therefore  use  the  FFT's  to  invert  (9)  to  obtain  the  values  of  the  A^s.  We 
may  then  represent  the  spatial  derivative  <£/  3K  as  a  Chefcyshev  series  given 
by: 


as 

ax 


N 

-  E  A 

n 


n-o 


(1) 


T„w 


(10) 


Because  of  properties  of  Chetyshev  polynomials  (equation  4b)  we  may  relate 
the  spectral  coefficients  of  the  spatial  derivative,  ,to  the  known 

spectral  coefficients  of  E,  namely  An,  by  the  following  recurrence  relation. 


(1) 

n 


2 

C 


N 

E 

p^n+1 

p-Ha^odd 


(ID 


Since  the  An's  are  known  at  the  current  time  step  tQ  (not  necessarily  zero) 
from  Equation  9,  the  A^^'s  are  obtained  from  the  recurrence  relation, 
Equation  11.  The  sunmation  in  Equation  10  is  performed  using  the  FFT. 
Therefore  it  remains  only  to  calculate  the  temporal  derivative  -g-  in  (7). 
For  the  results  presented  herein  the  Adams-Bashforth  algorithm  was  used  to 
advance  the  solution  to  t  +  £  •  (The  modified  Euler  predictor  corrector 
scheme  was  also  investigated.  However,  it  did  not  yield  better  results  and 
took  more  computer  time  to  implement.)  This  process  is  then  cyclically 
repeated  to  march  the  solution  in  time  (physical  or  coup u tat ional)  •  The 
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Adams -Bashforth  algorithm  is  given  by: 

r i  *  fgf‘  02) 


where  superscripts  denote  the  value  of  time  at  which  each  term  is  evaluated. 


3.  ONE  DIMENSIONAL  PROPAGATING  SHOCK  WAVE  RESULTS 

Two  types  of  artificial  viscosity  schemes  were  tried;  a  second  order 
scheme  given  by 

D  4  “  ~v(U  ix,  -  2U  4  +  u  .  .  )  (13) 

n,i  v  n,  i+1  n,i  n,i-l  J 


and  a  fourth  order  scheme  given  by 


n,i 


-w  to  ,  +  u  ,  -4  tu  ,  +  U  ,]  (14) 

1  n,  1+2  n,  1-2  n,i+l  n,  1-1 


+6 


Vi> 


where  Dn  ^  is  the  magnitude  of  the  dissipation  for  the  n^  conservative  flow 
variable  at  the  i spatial  point.  In  both  cases,  y  is  the  magnitude  of  the 
artificial  viscosity. 

Three  types  of  shock  tube  flews  were  considered:  (a)  supersonic 
inflow/outflow,  (b)  subsonic  inflow/outflow  and  (c)  supersonic 
inflow/subsonic  outflow.  They  represent  the  entire  range  of  shodt  tube 
problems  and  will  be  discussed  below.  The  time  step  size  used  throughout 
was  one  half  the  maximum  based  on  stability  considerations  (effectively  a 
Courant  nrnber  of  0.5).  The  resulting  time  step  size  values  are 
(a)  .502xl0-4,  (b)  .99£bcl0-4  and  (c)  .574xl0"4. 
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The  conditions  for  the  supersonic  inflow/ou  tflow  case  were  a  free 
stream  Mach  number  of  1.5  and  a  shock  Mach  number  of  3.5  (with  respect  to 
ground  fixed  coordinates).  One  hundred  twenty  eight  Chetyshev  polynomial 
terms  were  used  to  represent  the  flow.  All  results  utilized  the  low  pass 
spectral  filter.  In  all  cases  the  initial  shock  position  (t*tQ*0) 
was  at  x  *  -1.0  (i.e.  at  the  left  hand  side  conputational  boundary). 

The  second  order  artificial  viscosity  scheme  was  used  first  for  the 
above  problem.  Typical  results  are  shown  in  Figure  2  where  the  static 
pressure  distribution  (non-dime nsionalized  by  the  free  stream  va  i  p^)  is 

shown.  ?2  represents  the  post  shock  static  pressure.  The  anal:  ’  shock 

wave  position  at  t».1505  is  shown  for  comparison.  Clearly  the  ?  *  wave  is 

unacceptably  smeared.  For  this  reason  the  second  order  smoothii  eme  was 

abandoned. 

Results  for  this  case  with  the  fourth  order  smoothing  are  shown  in 
Figures  3  and  4.  The  figures  show  the  calculated  shock  solution  at  times  of 
0.05  and  0.15  respectively.  The  analytic  shock  position  at  the  respective 
times  is  shown  for  comparison.  As  can  be  seen,  the  confuted  shock  position 
is  in  excellent  agreement  with  the  analytic  solution.  Further,  the  correct 
pre  and  post  shock  pressures  are  maintained.  One  can  see  the  effect  of  grid 
resolution  by  conparing  Figures  3  and  4.  As  previously  mentioned  in  Section 
2  points  are  clustered  about  x  ■  i  1  with  the  coarsest  grid  spacing 
occurring  at  x-0.  The  shock  wave  is  in  a  region  of  high  point  resolution  in 
Figure  3  and  nearly  at  the  most  coarse  grid  resolution  in  Figure  4.  The 
apparent  skewness  of  the  calculated  shock  front  in  Figure  4  is  not  due  to 
overly  large  dissipation.  It  is  instead  due  to  the  coarse  grid  spacing. 
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The  shock  cannot  of  course  be  resolved  to  within  a  single  grid  spacing.  All 
flow  properties  were  held  fixed  at  both  the  supersonic  inflow  and  outflow 
boundaries.  At  the  supersonic  outflow  boundary  it  was  necessary  to  apply 
the  second  order  artificial  viscosity  locally  in  order  to  remove 
oscillations  emanating  from  this  boundary.  Without  this  localized  second 
order  smoothing,  the  solution  went  catastrophically  unstable  at  the  outflow 
boundary . 

The  second  shock  tube  problem  considered  involved  supersonic  inflow  and 
subsonic  outflow  (see  Figures  5  and  6).  The  free  stream  Mach  number  was 
0.  845  with  the  shock  Mach  number  2.949  with  respect  to  the  ground.  Again 
the  shock  is  maintained  as  a  sharp  discontinuity  propagating  at  the  correct 
velocity.  As  before,  the  supersonic  inflow  boundary  conditions  are  all  flow 
variables  held  fixed.  However,  at  the  subsonic  outflow  boundary  one 
physical  flow  variable  was  specified  with  the  remaining  ones  conputed  from 
the  characteristic  values  of  the  flow  (as  in  Reference  2). 

The  last  shock  tube  problem  considered  had  both  subsonic  inflcv  and 
subsonic  outflow.  The  free  stream  mach  nunber  was  0.5  while  the  shock  wave 
Mach  number  was  1.8  with  respect  to  the  ground.  Results  for  the  pressure 
distribution  at  two  different  times  are  shown  In  Figures  7  and  8.  As  in 
previous  cases,  the  shock  position  and  shape  are  in  excellent  agreement  with 
the  analytical  values.  The  boundary  conditions  used  were  to  hold  all  flow 
variables  fixed  at  the  subsonic  inflow  boundary  »nd  (as  in  the  previews 
case)  to  hold  one  flow  variable  fixed  at  the  subsonic  outflow  boundary  while 
conputing  the  remaining  ones  from  the  characteristics. 
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4.  TWO  DIMENSIONAL  SUPERSONIC  WEDGE  FLOW  RESULTS 


Two  cases  were  considered,  a  ten  degree  half  angle  wedge  at  free  stream 
Mach  numbers  of  1.5  and  3.0.  The  computational  grid  was  dimensioned 
33  x  33.  The  resulting  grid  lines  are  plotted  in  Figure  9.  Figure  10  shows 
the  allignment  of  the  computational  boundary  in  physical  space.  Now,  since 

X1  *  X  <  X2 

(15) 


we  must  transform  to  (  £,  n  )  space  to  obtain  |c|<  1,  ln!<  1  , which  is  required 
of  the  collocation  points.  This  transformation  is  given  by: 


5  - 


n  - 


2x-(x1-hc2) 

2y-y 

J  7  max 
y 

max 


(16) 


No  attempt  was  made  to  use  the  optimim  time  step  size,  given  by  Reference  5. 


(*) 


max 


8.0 


N  u+c 


max 


(17) 


In  fact,  for  all  calculations  presented  herein,  an  effective  Courant  number 
of  0.5  was  used.  (That  is,  the  numerator  of  (17)  was  replaced  by  4.0.)  For 
purposes  of  comparison,  the  wedge  surface  pressure  and  density  distributions 
as  well  as  computer  generated  contour  plots  of  the  shock  wave  position  and 


shape  are  used.  Fourth  order  dissipation  was  used  throughout  in  both  the 
x  and  y  directions*  Second  order  dissipation  was  used  in  the  neighborhood 
of  x  «  (in  the  x-direction  only)  to  eliminate  oscillations  emanating  from 
the  supersonic  outflow  boundary*  T^«  flow  Internal  to  the  conputational 
boundaries  was  initialized  to  free  stream  values*  Along  region  BCDE  of  the 
cooputational  boundary  the  flow  was  held  fixed  at  free  stream  values.  Along 
region  AB  and  FE  it  was  held  fixed  at  wedge  flow  properties*  Finally  at  the 
wedge  surface,  region  AF,  surface  tangency  was  imposed  after  each  time  step* 

Results  for  the  Mach  1.5  case  are  shown  in  Figures  11  through  14*  The 
time  step  size  was  .125x10  •  Figures  11  and  12  show  contour  plots  of  the 

pressure  and  density  fields  respectively*  The  analytic  shock  position  is 
also  shown  as  a  solid  line  for  .comparison.  The  shock  position  and 
orientation  are  predicted  exactly  by  the  pseud  os  pectral  solution*  The  wide 
band  or  thickness  of  the  computed  shock  is  due  to  the  very  coarse  grid 
resolution  used  in  the  2D  runs,  namely  33x33.  In  terras  of  grid  intervals 
the  shock  shown  In  Figures  11  and  12  lies  over  only  two  to  three  grid 
intervals*  Increasing  the  grid  resolution  will  reduce  the  thickness  of  the 
contoured  shock  wave. 

The  increase  in  shock  thickness  which  appears  in  Figures  11  and  12  in 
the  neighborhood  of  the  right  hand  side  cooputational  boundary  is  due  to  the 
localized  second  order  artificial  viscosity  scheme  used  in  that  region 
(supersonic  outflow).  With  a  33  point  resolution  along  the  x-axis  only 
several  points  are  needed  to  physically  extend  well  into  the  interior  of  the 
cooputational  area*  With  a  more  realistic  grid  resolution,  say  128  points, 
the  maximum  extent  would  be  rediced  to  only  x  *  0.99  and  no  effect  would  be 
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present  in  the  shock  plots*  Surface  pressure  and  density  distributions  are 
shown  in  Figures  13  and  14*  In  both  cases,  agreement  between  the  coaputed 
result  and  the  analytic  result  (represented  by  the  dotted  line)  is 
excellent.  The  minor  overshoots  and  undershoots  that  appear  in  both  plots 
represent  differences  of  less  that  1.5%  from  the  analytic  values. 

Similar  plots  are  shown  in  Figures  15  through  18  for  the  Mach  3.0 
case.  The  time  step  size  is  .  785x10  .  Agreement  is  excellent  both  in  the  shock 
shape  and  location  and  in  the  surface  pressure  and  density  distributions. 

5.  CONCLUSIONS 

(1)  Pseud ospectral  solution  techniques  can  treat  inviscid  l-D  and  2-D 
flows  with  shock  waves  quite  accurately  when  a  low  pass  spectral  filter  is 
used  in  conjunction  with  a  fourth  order  artificial  viscosity  scheme  (applied 
to  the  physical  variables).  All  shocks  are  maintained  as  discontinuities 
with  only  minor  pre  and  post  cursor  oscillations. 

(2)  For  the  2D  flow  problem  considered  here  (as  well  as  the  ID 
supersonic  outflow  problem),  a  localized  second  order  artificial  viscosity 
scheme  oust  be  applied  in  the  neighborhood  of  the  supersonic  outflow 
boundary  to  damp  out  oscillations  that  arise  at  the  boundary  and  keep  the 
solution  stable.  Without  it,  the  solution  always  goes  catastrophically 
unstable  at  this  boundary. 
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Fig.  1  —  Pseudospectral  calculation  flowchart 
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COMPUTATIONAL  BOUNDARY 
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Fig.  13  —  Comparison  of  analytic  and  calculated  wedge  surface  pressure  distributions 


Pig.  14  —  Comparison  of  analytic  and  calculated  wedge  surface  density  distributions 


Fig.  15  —  Supersonic  wedge  flow  computed  shock  wave 
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